#load pakages

if(!require(pacman)){install.packages("pacman")}
pacman::p_load(tidyverse, haven, ggpubr, rstatix, glue, estimatr, ggplot2, MetBrewer)


#load data
data <- read_dta("data/hgopy_pilot_public.dta")
#create the new price variable with discounts grouped and counselling style as factors
data <- data %>% mutate(data, 
                        pdx = ifelse(price_larc_cat==3,'Full','Discounted'),
                        pdn = ifelse(price_larc_cat==3,0,1), 
                        pd = factor(pdx),
                        pd = fct_reorder(pd,pdn),
                        csx = ifelse(view==0,'MiM',
                                     ifelse(view==1,'IDM',
                                            ifelse(view==2,'SDM','None'))),
                        cs = factor(csx),
                        cs = fct_reorder(cs,view),
                        price_larc = fct_relevel(as.character(price_larc), "0", "150", "1000", "2000", "4000"),
                        price_larc = case_when(
                          price_larc %in% c("4000") ~ "Full",
                          price_larc %in% c("2000", "1000","150") ~ "Discounted",
                          price_larc == "0" ~ "Free"
                        ),
                        price_larc = fct_relevel(price_larc, "Free", "Discounted", "Full")
)


###Plot prep 

# define outcome for convenience
data = data %>% mutate(y=adopted_larc)

#tests for the within group comparisons, Only MiM
tests_wncs_mim <- data %>%
  select(cs,price_larc,y) %>%
  filter(cs=='MiM') %>% 
  group_by(cs) %>% 
  t_test(y ~ price_larc, detailed = TRUE)
#gotta add the y position...
tests_wncs_mim <- tests_wncs_mim %>% 
  add_xy_position(x = "cs" , dodge = 0.8) # , fun = "mean_se"
#format the p values and estimates nicely 
tests_wncs_mim <- tests_wncs_mim %>% 
  mutate(prounded = round(p, digits=3),
         est = round(estimate, digits=3))
print(tests_wncs_mim)

#tests for the within group comparisons, excluding MiM
tests_wncs_idm <- data %>%
  select(cs,price_larc,y) %>%
  filter(cs=='IDM') %>% 
  group_by(cs) %>% 
  t_test(y ~ price_larc, detailed = TRUE)
#gotta add the y position...
tests_wncs_idm <- tests_wncs_idm %>% 
  add_xy_position(x = "cs" , dodge = 0.8) # , fun = "mean_se"
#format the p values nicely 
tests_wncs_idm <- tests_wncs_idm %>% 
  mutate(prounded = round(p, digits=3),
         est = round(estimate, digits=3))
print(tests_wncs_idm)

#tests for the within group comparisons, in SDM
tests_wncs_sdm <- data %>%
  select(cs,price_larc,y) %>%
  filter(cs=='SDM') %>% 
  group_by(cs) %>% 
  t_test(y ~ price_larc, detailed = TRUE)
#gotta add the y position...
tests_wncs_sdm <- tests_wncs_sdm %>% 
  add_xy_position(x = "cs" , dodge = 0.8) # , fun = "mean_se"
#format the p values nicely 
tests_wncs_sdm <- tests_wncs_sdm %>% 
  mutate(prounded = round(p, digits=3),
         est = round(estimate, digits=3))
print(tests_wncs_sdm)


###Plot - definition 1, without LAM 
ggbarplot(data, x = "cs" , y = "y", 
          fill = "price_larc", 
          add = c("mean_ci"),
          position = position_dodge(0.8),
          xlab = "Counselling style" , 
          ylab = "", # "Fraction of clients",
          title = "Fraction of clients who adopted a LARC", 
          label = TRUE,
          lab.nb.digits = 3,
          lab.hjust = -.35,
          lab.size = 3,
          ylim = c(0,1.3)) +
  #add in the p-values and differences, all default set at 1 so have to manually lower them down
  stat_pvalue_manual(tests_wncs_idm, label = "d = {est}; p = {prounded}",
                     tip.length = 0.005,
                     step.increase = 0.15,
                     bracket.nudge.y = -.15,
                     label.size = 3) +
  stat_pvalue_manual(tests_wncs_sdm, label = "d = {est}; p = {prounded}",
                     tip.length = 0.005,
                     step.increase = 0.15,
                     bracket.nudge.y = -.15,
                     label.size = 3) +
  stat_pvalue_manual(tests_wncs_mim, label = "d = {est}; p = {prounded}",
                     tip.length = 0.005,
                     step.increase = 0.15,
                     bracket.nudge.y = -.15,
                     label.size = 3) +
  #make it look nice
  theme_bw() +
  scale_fill_manual(values =met.brewer("Egypt", 3)) + #c( "#43b284", "#0f7ba2", "#dd5129") 
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=12),
        legend.text = element_text(size=12),
        legend.position= "bottom",
        legend.title = element_text(size=12),
        panel.grid.minor = element_blank(), 
        panel.grid.major = element_blank(),
        plot.caption = element_text(hjust=0,size=10)) +
  guides(fill=guide_legend(title="LARC prices:")) 

ggsave("output/figures/FigS1-view-prices-disaggregated-cat.png", dpi = 320, scale = 0.8, height = 16, width = 24, units = "cm")
